#########################
##### Punish or Tolerate? 
##### Replication File
#########################

## Packages 
library(tidyverse)
library(lmtest)
library(stargazer)
library(sandwich)
library(AER)

rm(list=ls())

setwd("~/Dropbox/sexual violence/data/")

## Open Dataset 
master <- read.csv("master.csv")


###########################
############## Descriptive Statistics: Table 2 
###########################

descr <- master %>% dplyr::select(sv_mean, independent, ccs_dynamic, active_count, 
                                  log_trade, polityiv, gwf_military, v2clgencl, loggdp, 
                                  logexc, logpop, logbrd, reb_n, duration, pcount)
descr <- as.data.frame(descr)
stargazer(descr, title = "Descriptive Statistics", digits = 1,
          covariate.labels = c("Avg. Sexual Violence", "Independent","Civilian Control",
                               "Number of Militias",
                               "Trade (% GDP)", "Democracy Score", "Military Regime", "Gender Equality", 
                               "GDP per capita", "Excluded Groups", "Population",
                               "Battle Deaths",  "Number of Rebels", "Duration",   "Number of Security Forces" ))

###########################
############## Descriptive Values of Sexual Violence by Civilian Control: Figure 1
###########################

summary(master$ccs_dynamic)
master$ccs_cat <- ifelse(master$ccs_dynamic <= -0.79893, "Low", 
                         ifelse(master$ccs_dynamic >= 0.06142, "High", "Mod"))

table(master$ccs_cat)
submaster <- master %>% group_by(ccs_cat) %>% summarize(mean_sgbv = mean(sv_mean, na.rm = T))
submaster <- submaster %>% filter(!(is.na(ccs_cat)|ccs_cat == "Mod"))

submaster$ccs_cat <- factor(submaster$ccs_cat, levels = c("Low", "High"))

pdf("descr_fig.pdf",height=6.1,width=6.6)
ggplot(data = submaster, aes(x = ccs_cat, y = mean_sgbv)) + 
  geom_bar(stat = "identity", aes(fill = ccs_cat)) + 
  scale_fill_grey(start = 0.5, end = .7) + 
  theme_bw() +
  xlab("Civilian Control") + 
  ylab("Average Level of Sexual Violence") +
  theme(legend.position = "none")
dev.off()


###########################
############## OLS Models: Table 3
###########################

mod1 <- lm(sv_mean ~ independent + log_trade + polityiv + gwf_military +  v2clgencl + loggdp  
           + logexc + logpop + logbrd  + reb_n + duration + pcount + as.factor(year), data = master)
m1s <- coeftest(mod1, vcov = vcovHC(mod1, type="HC1"))

mod2 <- lm(sv_mean ~ ccs_dynamic   +  log_trade +polityiv + gwf_military  +  v2clgencl + loggdp 
           + logexc + logpop + logbrd  + reb_n + duration +pcount + as.factor(year), data = master)
m2s <- coeftest(mod2, vcov = vcovHC(mod2, type="HC1"))

mod3 <- lm(sv_mean ~ active_count  +  log_trade + polityiv  + gwf_military  +  v2clgencl + loggdp  
           + logexc + logpop +  logbrd  + reb_n + duration + pcount + as.factor(year), data = master)
m3s <- coeftest(mod3, vcov = vcovHC(mod3, type="HC1"))


stargazer(mod1,mod2,mod3, se=list(m1s[,2],m2s[,2],m3s[,2]),type = "text")
stargazer(mod1,mod2,mod3,se=list(m1s[,2],m2s[,2],m3s[,2]),
          covariate.labels = c( "Independent","Civilian Control",
                                "Number of PGMs",
                                "Trade (% GDP)", "Democracy Score", "Military Regime", "Gender Equality", 
                                "GDP per capita", "Excluded Groups", "Population",
                                "Battle Deaths",  "Number of Rebels", "Duration", "Number of Security Forces"),
          dep.var.labels = "Average Level of Wartime Sexual Violence")


###########################
############## Predicted Values: Figure 2 
###########################

getmode <- function(v) {
  uniqv <- unique(v)
  uniqv[which.max(tabulate(match(v, uniqv)))]
} #function set up to get mode values for dichotomous variables 

sim_data <- data.frame(
  ccs_dynamic=seq(from=min(master$ccs_dynamic, na.rm = T),to=max(master$ccs_dynamic,na.rm=T),length.out=100),
  log_trade = mean(master$log_trade, na.rm = T), 
  polityiv = mean(master$polityiv, na.rm = T), 
  gwf_military= getmode(master$gwf_military),
  loggdp =  mean(master$loggdp,na.rm=T), 
  v2clgencl = mean(master$v2clgencl,na.rm=T), 
  logexc = mean(master$logexc,na.rm=T), 
  logpop = mean(master$logpop, na.rm = T), 
  logbrd = mean(master$logbrd, na.rm = T),
  reb_n = mean(master$reb_n, na.rm = T), 
  duration = mean(master$duration, na.rm = T),
  pcount = mean(master$pcount, na.rm = T),
  year = 2002)
OutHats<-predict(mod2, newdata=sim_data, se.fit = TRUE)
OutHatsUB<-OutHats$fit+(1.96*OutHats$se.fit)
OutHatsLB<-OutHats$fit-(1.96*OutHats$se.fit)
OutHats<-cbind(as.data.frame(OutHats),OutHatsUB,OutHatsLB)
both<-cbind(sim_data,OutHats)

pdf("pred_fig2.pdf",height=6.1,width=6.6)
plot(NULL,
     xlim = c(min(master$ccs_dynamic, na.rm= T),max(master$ccs_dynamic, na.rm = T)),
     ylim = c(0, 3),
     axes = F, xlab = NA, ylab = NA) 
polygon(c(both$ccs_dynamic,rev(both$ccs_dynamic)),c(both$OutHatsUB,rev(both$OutHatsLB)),col="gray95",border="gray95")
grid()
lines(both$ccs_dynamic,both$fit,col="black",lwd=1.75)
points(jitter(master$ccs_dynamic),rep(0,length(master$ccs_dynamic)), pch="|",col=alpha("#000000",.7))
box()
axis(1)
axis(2)
mtext(side=1,"Civilian Control",line=2)
mtext(side=2,"Predicted Level of Sexual Violence",line=2.5)
title("")
dev.off()

###########################
############## Tobit Models: Table 4
###########################

#VCov HC does not run in tobit models, instead we use clustered robust errors 
#https://stackoverflow.com/questions/66412110/sandwich-mlogit-error-in-ef-x-non-conformable-arrays-when-using-vcovhc

t1 <- tobit(sv_mean ~ independent  + log_trade +polityiv + gwf_military  +  v2clgencl + loggdp 
            + logexc + logpop + logbrd  + reb_n + duration + pcount+ as.factor(year),   
            data = master)
t1cl <- coeftest(t1, vcov = vcovCL)

t2 <- tobit(sv_mean ~ ccs_dynamic   + log_trade +polityiv + gwf_military  +  v2clgencl + loggdp 
            + logexc + logpop + logbrd  + reb_n + duration + pcount + as.factor(year),   
            data = master)
t2cl <- coeftest(t2, vcov = vcovCL)

t3 <- tobit(sv_mean ~ active_count   + log_trade +polityiv + gwf_military  +  v2clgencl + loggdp 
            + logexc + logpop + logbrd  + reb_n + duration + pcount + as.factor(year),   
            data = master)
t3cl <- coeftest(t3, vcov = vcovCL)

stargazer(t1,t2,t3,se=list(t1cl[,2],t2cl[,2],t3cl[,2]), type = "text")

stargazer(t1,t2,t3,se=list(t1cl[,2],t2cl[,2],t3cl[,2]), 
          covariate.labels = c( "Independent", "Civilian Control",
                                "Number of PGMs", 
                                "Trade (% GDP)", "Democracy Score", "Military Regime", "Gender Equality", 
                                "GDP per capita", "Excluded Groups", "Population",
                                "Battle Deaths",  "Number of Rebels", "Duration",   "Number of Security Forces"),
          dep.var.labels = "Average Level of Wartime Sexual Violence")

